########################################################
# Replication file for Ciudad Autonoma de Buenos Aires #
####################### Turnout ########################
########################################################
rm(list=ls(all=TRUE))
#setwd("")
library(eiPack)
library(coda)
library(foreign)

data=read.csv(file.path(getwd(),"data","CABA_data.csv"))


# Define row and column marginals of the transition matrix. 
UP=data$paso_9002
OP=data$votantes-data$paso_9002
UG=data$gral_9002
OG=data$votantes-data$gral_9002

data1=data.frame(UP,OP,UG,OG)


# Tunning the EI algorithm
tune.nocov=tuneMD(cbind(UG,OG)~cbind(UP,OP),data=data1,ntunes=1,
	     totaldraws=200000)
# MCMC
out.nocov=ei.MD.bayes(cbind(UG,OG)~cbind(UP,OP),covariate=NULL,
          data=data1,tune.list=tune.nocov,verbose=0,
	    ret.beta="r",thin=80,burnin=100000)  

# Obtain transition matrix for each mesa: 
source(file.path(getwd(),"Codes","getbetas.R"))
betaUPUG=getbetas(out.nocov,"UP","UG")
betaUPOG=getbetas(out.nocov,"OP","UG")
betaOPUG=getbetas(out.nocov,"OP","OG")
betaOPOG=getbetas(out.nocov,"UP","OG")

betaout=data.frame(betaUPUG,betaUPOG,betaOPUG,betaOPOG,data$votantes,
                   data$comuna,data$mesa)
write.table(betaout,file=file.path(getwd(),"Output","TurnoutCABA_Mesa.txt"),col.names=NA)



# Obtain aggregate transition matrix. 
source(file.path(getwd(),"Codes","getbetasagr.R"))
betaUPUG=getbetasagr(out.nocov,"UP","UG")
betaOPUG=getbetasagr(out.nocov,"OP","UG")
betaOPOG=getbetasagr(out.nocov,"OP","OG")
betaUPOG=getbetasagr(out.nocov,"UP","OG")
betaoutagr=data.frame(betaUPUG,betaUPOG,betaOPUG,betaOPOG,data$votantes,
                   data$comuna,data$mesa)

betaoutagr=as.vector(betaoutagr[1,1:16])
betaoutagr=matrix(betaoutagr,4)
betaoutagr_coefs=t(matrix(betaoutagr[1,],2))
betaoutagr_sd=t(matrix(betaoutagr[2,],2))
write.table(betaoutagr_coefs,file=file.path(getwd(),"Output","Turnout_CABAAggregate_Coefs.txt"),col.names=NA)
write.table(betaoutagr_sd,file=file.path(getwd(),"Output","Turnout_CABAAggregate_SD.txt"),col.names=NA)



# Geweke Convergence Statistics
geweke=geweke.diag(out.nocov$draws$Beta)
gewekestats=geweke[[1]]
gewekedf=data.frame(gewekestats)
#write.table(gewekedf,file="gewekelist.txt",col.names=NA)


jpeg(file.path(getwd(),"Output","Turnout_CABAHistogramGeweke.jpeg"))
hist(gewekestats)
dev.off()

print("Max Geweke")
print(max(gewekestats))
print("Min Geweke")
print(min(gewekestats))
print("Percentage Above 2.35")
print(length(gewekestats[abs(gewekestats)>2.35])/length(gewekestats))
print("Percentage Above 1.96")
print(length(gewekestats[abs(gewekestats)>1.96])/length(gewekestats))
print("Percentage Above 1.5")
print(length(gewekestats[abs(gewekestats)>1.5])/length(gewekestats))
print("Percentage Above 1.5")
print(length(gewekestats[abs(gewekestats)>1])/length(gewekestats))
print("Percentage Above 0.5")
print(length(gewekestats[abs(gewekestats)>0.5])/length(gewekestats))
